%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PNGraphene(N,L) is a function that computes the charge and potential
% distribution of a Graphene p-n junction. 
% Input: N, L are the Grid Size and the Periodic Length, respectively.
% Output: V,R are the Potential and Charge distribution, in real space.
% Note: All calculations are done in Fourier Space
%


function [V,R,Err] = PNGraphene(f,mf,N,L)
%% Function Setup
global a b c q;         % Constant Coefficients Declared as Global 
% Fourier Grid Setup
if mod(N,2) == 0        % Making the size N an odd number to make the k 
    N = N + 1;          % grid symmetric about the origin
end
n = -(N-1)/2:(N-1)/2;
oddk_index = 2:2:N-1;     % index of the odd k points
evenk_index = 1:2:N;  % index of the even k points
nc = (N+1)/2;
h = L/(N-1);            % grid size in real space (in Bohr Radius)

% Physical Constant Setup (In Hartree Atomic Unit)
eps = 1;            % Relative permittivity
eps_0 = 1/(4*pi);   % Vacuum permittivity 
%A = 0.6551/(0.5291772109217)^2; % Area per atom

vf = 1E6/(2.187691263373E6); 
e = 1; % Electronic Charge
hbar = 1;
kn = (2*pi/L)*n;


% Equation Coefficients
a = (e*f/eps);
b = ((e^3))/(eps*pi*(hbar^2)*(vf^2));
c = 1/(2*eps*eps_0);
%mf = 1E-2;    % Potential Mixing Factor

% Square Wave in k space (for symmetry)
q = zeros(1,N);
q(oddk_index) = (1/pi)*(2/1i)*(1./n(oddk_index));
q(nc) = 0;
%figure; plot(N*real(ifft(ifftshift(conv(q,conv(q,q,'same'),'same')))));


% Initial Setup for v0 in k space
Vmax = sqrt(pi*f*(hbar^2)*(vf^2)/(e^2));
v0 = Vmax*q;
% figure; plot(real(ifft(ifftshift(v0))));
% Error 
err = Inf;iter=0;
Err=[];
%%%%%%%%%%%%%%%%%%%%%%%% End of Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Function Body
while err > 1E-6
    
    % v-> grDOS -> s
    s = grDOS(v0);
    % s -> Gauss' Theorem -> v
    v = c*s./abs(kn);
    v(nc) = 0;
    
    % Recording the error
    err = sum(abs(v-v0))*h;
    Err = [Err,err];
    % mixing
    v0 = (1-mf)*v0 + mf*v;    
    display(sprintf('err(%d)=%.6e',iter,err));
    iter=iter+1;
end

V = N*fftshift(ifft(ifftshift(v)));
R = N*fftshift(ifft(ifftshift(s)));








%% My Functions
%%%%%%%%%% Function grDOS %%%%%%%%%%
    function Sk = grDOS(Vk)
        Vk1 = conv(Vk, Vk,'same');
        Vk2 = conv(q, Vk1,'same');
        Sk = a*q - b*Vk2;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%% Function MyDraw %%%%%%%%%%
    function MyDraw(Vx, Sx, fig_handle)
        
        figure(fig_handle);
        subplot(2,1,1);  
        plot(real(Vx)); hold on; 
        plot(ones(1,N)*sqrt(pi*f*(hbar^2)*(vf^2)/(A*e^2)),'--k');
        plot(-ones(1,N)*sqrt(pi*f*(hbar^2)*(vf^2)/(A*e^2)),'--k');
        ylabel('V(x)');
        hold off;
        
        subplot(2,1,2); 
        plot(real(Sx));ylabel('\sigma(x)');
        
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        

end